Radiative Falloff in Neutron Star Spacetimes 
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We systematically study late-time tails of scalar waves propagating in neutron star spacetimes. 
We consider uniform density neutron stars, for which the background spacetime is analytic and 
the compaction of the star can be varied continously between the Newtonian limit 2M/R <C 1 and 
the relativistic Buchdahl limit 2M/R — 8/9. We study the reflection of a finite wave packet off 
neutron stars of different compactions 2M/R and find that a Newtonian, an intermediate, and a 
highly relativistic regime can be clearly distinguished. In the highly relativistic regime, the reflected 
signal is dominated by quasi-periodic peaks, which originate from the wave packet bouncing back 
and forth between the center of the star and the maximum of the background curvature potential 
at R ~ 3M. Between these peaks, the field decays according to a power-law. In the Buchdahl 
limit 2M/R — » 8/9 the light travel time between the center and the maximum or the curvature 
potential grows without bound, so that the first peak arrives only at infinitely late time. The modes 
of neutron stars can therefore no longer be excited in the ultra-relativistic limit, and it is in this 
sense that the late-time radiative decay from neutron stars looses all its features and gives rise to 
power-law tails reminiscent of Schwarzschild black holes. 
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I. INTRODUCTION 

The late-time radiative falloff in a black hole space- 
time was first studied with the goal of understanding 
how, in the absence of rotation, a spherically symmetric 
Schwarzschild black hole emerges in a general aspheri- 
cal collapse. The first numerical results by de la Cruz, 
Chase and Israel jl| were followed up by the perturba- 
tive analytic study by Price who showed that scalar, 
electromagnetic and gravitational radiation all decay as 
an inverse power of time. This latc-timc behavior is com- 
monly refered to as "power-law tails" . 

The analysis of Gundlach, Price and Pullin H showed 
that power-law tails arise from backscattering off the 
weak gravitational potential in the far zone (c/. One 
immediate implication, namely that the late time behav- 
ior should be different in black hole spacetimes which are 
not asymptotically flat, was confirmed by Brady, Cham- 



bers, Krivan and Laguna 
Laarakkers and Poisson 



and later Brady, Chambers, 
These authors studied ra- 
diative falloff in Schwarzschild-de Sitter spacetimes, and 
found that for nonspherical (and minimally coupled) per- 
turbations, the power-law decay at early times changed 
into an exponential decay after a cosmological horizon 
timescale . Another implication is that power-law tails 
should be present even in the absence of event horizons, 
e.g. for perturbations of neutron stars. 

Perturbations of nonrotating, spherically symmetric 
relativistic neutron stars have been studied by many au- 
thors, mostly with the goal of finding their oscillation 
modes [pHlo]- Typically, these modes are determined by 
assuming a time-dependence exp(zert) and then studying 
the time-independent problem. The perturbation equa- 
tions then reduce to an eigenvalue problem, which can 
be solved for the complex frequencies a. In addition to 



families of modes corresponding to those of Newtonian 
stars, relativistic stars also have modes which are remi- 
niscent of the quasi-normal modes of black holes and can 
be associated with the dynamical degrees of freedom of 
the gravitational field. Because of their association with 
gravitational waves, these modes are usually refered to 
as w-modes. Several distinct families of w-modes, both 
in axial and polar parity, have been found in numeri- 
cal calculations (see the recent review by Kokkotas and 
Schmidt, |ft7j] ): curvature modes exist for all compactions 
2M I R, where M is the stellar mass and R the circum- 
ferential radius, trapped modes only arise in sufficiently 
compact stars with 2M/R > 2/3, when radiation can be 
"trapped" inside the maximum of the background cur- 
vature potential at R ~ 3M (see below), and interface 
modes are associated with a discontinuity of the back- 
ground curvature at the surface of the star. 

Similar to black holes, the response of a neutron star 
to a general perturbations will be a superposition of fluid 
and gravitational wave modes and a power-law tail. Since 
all modes decay exponentially, the late-time falloff should 
be dominated by the power-law tail. Unlike the quasi- 
normal modes of a black hole, however, the frequency of 
the neutron star's w-modes depends on the compaction 
2M/R of the star. As the compaction approaches the 
ultrarelativistic Buchdahl limit 2M/R = 8/9 p|, the 
imaginary part of the frequency of at least some of these 
modes approches zero, implying that their damping time 
becomes infinite jt6|,[l7). A similar behavior had ear- 
lier been identified by Detweiler for the fluid modes (/- 
modes) of relativistic stars ||. This observation imme- 
diately raises the question how a neutron star in the ul- 
trarelativistic limit responds to a general perturbation, 
whether a power-law tail emerges, and whether its char- 
acteristic radiative fall-off approaches that of a black 
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hole. In this paper we seek to answer this question by 
dynamically evolving the time-dependent perturbation 
equations. 

So far, few authors have studied perturbations of 
relativistic stars by dynamically evolving the time-de- 
pendent perturbation equations, as opposed to assum- 
ing a exp(zcri) time-dependence and solving the time- 
independent equations. Andersson and Kokkotas Jigj] 
studied the excitation of w-modes by an impinging grav- 
itational wave, and discussed the prospect of detections 
by future gravitational- wave laser interferometers. Allen, 
Andersson, Kokkotas and Schutz |2(J analyzed the spec- 
trum of excited fluid and wave modes arising from several 
different sets of initial data. In some of these simulations, 
both the quasi-normal modes and subsequent power- 
law tails were observed (see also Fig. 5 in Jill). 
Ruoff generalized earlier results and included real- 
istic equations of state. 

In this paper, we dynamically evolve a scalar field in 
a neutron star background. We model the neutron star 
as a homogeneous, incompressible sphere, for which the 
background metric is known analytically. We derive an 
analytical generalized tortoise coordinate, which brings 
the equations into a particularly simple form. We then 
vary the compaction of the star between the Newtonian 
limit 2M/R < 1 and the Buchdahl limit 2M/R = 8/9. 
Interestingly, the generalized tortoise coordinate distance 
between the center and the surface of the star diverges 
in the Buchdahl limit, implying that the light travel time 
becomes infinite. 

We then study the reflection of finite wave packets off 
neutron stars of various compactions, and identify three 
distinct regimes: In the Newtonian regime 2M/R <C 1, 
we find a reflection of the initial pulse from the cen- 
ter of the star, and a subsequent power-law falloff. In 
an intermediate regime, so-called trapped modes of the 
neutron star are excited and decay exponentially (see 
also ]l9|-^l[|). After these modes have decayed to suf- 
ficiently small values, the late-time tail re-emerges as a 
power-law. In this paper we focus on the highly rela- 
tivistic regime, in which the wave packet bounces back 
and forth between the center of the star and the maxi- 
mum of the background curvature potential at R ~ 3M, 
giving rise to quasi-periodic peaks in the reflected signal. 
Between these peaks, the field again decays according to 
a power-law. In the Buchdahl limit 2M/R —> 8/9 the 
light travel time between the center and the maximum 
or the curvature potential grows without bound, so that 
the first peak arrives only at infinitely late time. The 
modes of neutron stars can therefore no longer be ex- 
cited in the ultra- relativistic limit, and it is in this sense 
that tails from neutron stars loose all their features and 
give rise to power-law tails. 

The paper is organized as follows. In Sec. II we intro- 
duce the basic equations, starting with a general spher- 
ically symmetric spacetime, then specializing to homo- 
geneous neutron stars, and describing our numerical im- 
plementation of the equations. We discuss our numerical 



results and present characteristic examples for each of 
three different compaction regimes in Sec. III. In Sec. IV 
we briefly summarize our findings. We also include an 
appendix containing the derivation a conserved energy 
integral of the motion, which we use as a numerical check. 



II. BASIC EQUATIONS 

A. Scalar Waves in Spherical Symmetry 

In spherical symmetry, the line element can be written 

ds 2 = -e 2>s> dt 2 + e 2A dr 2 + r 2 (d6 2 + sin 2 9d(j) 2 ). (1) 

The wave equation for a massless and minimally coupled 
scalar field <p, 
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n<p = <p>° 



can then be written 



e ip ttt H - (r e (pA 



(V=99 af} V, a ) fl =0, 



(2) 



(3) 

where we have decomposed the angular part of (p into the 
spherical harmonics (pit, r, 8, <j>) = (pit, r)Yg m (8, <j>). 

Introducing a "generalized tortoise coordinate" r* sat- 
isfying 



dr* 
dr 



„A-# 



and substituting 

ip = rip 

into eq. (||) yields 

- Ipjt + 1p,r.r, = Vl{r) if), 

where the effective potential Vi{r) is given by 



V e (r) = e 



V r 2 



+ (*-A),, 



(4) 

(5) 
(6) 

(7) 



B. Uniform Density Neutron Stars 

We now specialize to a constant density neutron star, 
for which 



3/ _2M\ 1/2 _1/ 2Mr 2 ^ 1/2 
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r < R 
r>R 



(8) 
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and 



yields the familiar Schwarzschild tortoise coordinate 
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r <R 
r > R. 
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Here R denotes the surface radius and M denotes the 
gravitational mass of the star |2^] . 

It is useful to introduce the nondimensional quantities 



r 

M ' 



r* = 



M 



~_ t 
~ M 



and to define the compaction of the star as 
C 



2M 
~R~' 



(10) 



(11) 



Note that C is limited by the well-known Buchdahl limit 
C < 8/9 @. With these definitions, the wave eq. (§) 
now becomes 



(12) 



where we have also used Ve(f) — M 2 Ve(r). In the follow- 
ing we will drop the tildes again, and it will be understood 
that all variables are nondimensional. This is equivalent 
to setting M = 1. 

Adopting the notation of Chandrasekhar and Fer- 
rari @, 



y = y/1 - CV/4, yi = VT~C, 
the generalized potential Ve(r) can be written 
( iyi-y (t(l+l), 



(13) 



V t {r) 



-(3yi - y) 



C 3 

f — (2y-3 yi 



r < 2/C 
r > 2/C. 



(14) 



In the exterior of the star r > 2/C the potential reduces 
to the Schwarzschild potential. Note also that the poten- 
tial is discontinous at the surface of the star (r = 2/C). 
The size of the discontinuity is 3C 3 (1 — C)/8, indepen- 
dent of I. In Figure (1) we show Vt{r) for different values 
of the compaction for £ = and £ = 1. 



r*=r + 21n(--l) 



(16) 



(where a constant of integration has been chosen appro- 
priately) . 

In the interior, r < 2/C, we have to integrate 



2dr 



V (3yi - y) ' 



With the substitution 



3j/i - y 

the integrand can be brought into the form 
4 f ds 



C 3 / 2 J (-y 2 s 2 + 6y 1 s-ly/ 2, 
where we have used 

2/2=8- 9C. 



(17) 



(18) 



(19) 



(20) 



The integration can now be carried out analytically and 
yields 



(C 3 y 2 )i/2 



112 



3yi - y 



3yi 



2yi 



3yi 



+ | + 21n(l-l), (21) 



where the dependence on r enters through y, and where 
we have chosen a constant of integration such that is 
continous across the surface of the star. As expected, 
we recover r*/r — > 1 in the Newtonian limit C — > 0. 
In the ultrarelativistic limit C — + 8/9, r* at the origin 
r = approaches negative infinity, r*(0) — + — oo. Since 
r* measures the coordinate time which elapses along a 
photon radial geodesic (eq. (12)), this implies that the 
light travel time between the surface and the center of 
the star grows without bound as the star's compaction 
reaches the Buchdahl limit. This property has important 
consequences for the scattering of waves from ultrarela- 
tivistic stars, as we will discuss in Sec. III. In Fig. 2, 
we show as a function of r for various different com- 
pactions C. 



C. The Generalized Tortoise Coordinate 

We now have to integrate eq. (Q) to find an expression 
for the generalized tortoise coordinate r*. For r > 2/C, 
the integration 

'-/■ M *-/i^f (15) 



D. Numerical Implementation 

We integrate eq. ( |l2|) numerically using the finite- 
differencing scheme described by Gundlach, Price and 
Pullin ||, which is based on the double null coordinates 

u = t — r*, v = t + . (22) 

In terms of these, eq. (p~2h becomes 
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FIG. 1. The potential inside and outside of the star for different values of the compaction C = 2M/R. Note the discontinuity 
at the stellar surface. 
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FIG. 2. The tortoise coordinate r* as a function of r for 
several different values of C. 



Denoting 



V e (r)ip. 



v + = v + A, 



(23) 



(24) 



where A is the grid spacing, we can write a finite differ- 
ence equation for eq. (E3J) as 



tp(u + ,v + )= ll-— ^(r)l (ip(u + ,v)+ip(u,v+)) 

-t/>(u,v). (25) 

Here r can be found from r* = (v — u)/2 by numerically 
inverting eq. (pl| ) . We impose a regularity boundary con- 
dition at the origin (tp = at r = 0) and an outgoing 
wave boundary condition 

tp{u+, v+) = ip{u + , v) (26) 

at a large radius r = r max in the asymptotically flat 
regime. 

We tested our code in various ways. For a vanish- 
ing potential Vi(r) (flat space and i = 0), the analytical 
scalar wave solution is reproduced exactly. For smooth 
potentials (as for a Schwarzschild black hole), the code 
produces second order accurate results (0(A 2 )). This 
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FIG. 3. Late-time behaviour for the £ — and 1 = 1 modes with C = 0.2. For such low values of the compaction we observe 
a power-law decay in a nearly Newtonian exterior geometry. 



second order convergence is broken by the discontinu- 
ity in the potential at the surface of the star, but the 
code is still convergent. We also convinced ourselves that 
by smoothing out the discontinuity, second order conver- 
gence is restored. Lastly, we checked that our evolution 
scheme conserves energy, which is discussed in more de- 
tail in Appendix A. 

Following ||, we choose as initial data ip(u, 0) = and 
ip(0,v) = exp {-(v - v c ) 2 /{2a 2 )) for v c = 30 and a = 3. 
However, the qualitative features of our findings should 
be independent of this particular choice. All results will 
be shown as a function of time t — ^(v+u) for an observer 
at fixed radius r* =50. 



III. NUMERICAL RESULTS 

We focus on the late-time radiative falloff of the scalar 
field, and find that three different regimes of neutron star 
compaction which can be clearly distinguished. We dis- 
cuss these different regimes below, and summarize our 
numerical results in Table 1. 

For small compactions we observe a power-law decay 



2M/R 


r.(0) 


r.(fl) 


a {I = 0) 


a {1=1) 


0.2 


0.704 


12.773 


-3.01 


-5.02 


0.83 


-16.799 


-0.762 


-3.01 


-5.02 


0.882 


-55.947 


-1.755 


-3.00 




0.888 


-163.039 


-1.889 


-3.01 




0.8888 


-525.380 


-1.907 


-3.01 





TABLE I. Summary of different runs. We tabulate the 
compaction C = 2M/R, the value of the tortoise coordinate 
at the origin r»(0) and the surface of the star r,(R), and the 
slope a of the late-time power law tail t~ a for I = and 1=1. 
For £ = with C = 0.888 and C = 0.8888, the slopes a refer 
to the power law decay before the arrival of the first of the 
periodic peaks. The accuracy of the slopes is limited by the 
finite integration time. 
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FIG. 4. Late-time behaviour for the I = and 1=1 modes with C = 0.83. Both modes are in the intermediate compaction 
regime where trapped quasi-normal modes are excited and are observed to decay exponentially. 



in a nearly Newtonian (exterior) Schwarzschild geome- 
try. The initial pulse is reflected at the origin and then 
propagates outward, leaving behind the power-law de- 
cay ip ~ t~ a from backscattering off the weak spacetime 
curvature at large radii. The slope a of this power-law is 
consistant with the value — (2£+3) predicted for backscat- 
tcr in any spherically symmetric, static, asymptotically 
flat spacetime || . As a typical example, we show results 
for C = 0.2 for both t = and I = 1 in Fig. 3. This 
behaviour is observed for compactions C up to about 0.8 
for the I — mode and 0.7 for the £ = 1 mode. 

At these high compactions, we find a smooth transi- 
tion into an intermediate regime. Once the star is suffi- 
ciently compact, a sufficiently deep potential well forms 
between the peak of the exterior Schwarzschild potential 
and the origin (see Fig. 1), and the incoming pulse can 
excite quasi-normal modes trapped in this potential well. 
These modes decay exponentially as their energy slowly 
"leaks" out to infinity. Once their amplitude has dropped 
to small enough values, the power-law tail emerges and 
dominates the late-time decay of the field. In Fig. 4, we 
demonstrate this behavior for C = 0.83. 

Note that due to the discontinuity of the potential, 



there is a small potential well just inside the surface of the 
star for all compactions C. For small compactions, how- 
ever, this well is too small to excite quasi-normal modes, 
and it is only when the surface of the star is close to the 
maximum of the exterior Schwarzschild potential that 
this potential well becomes sufficiently large (compare 
Fig. 1) to observe the exponential decay. 

A third kind of behaviour is observed in the strongly 
relativistic regime, when the compaction of the star is 
C > 0.86 for t = and C > 0.83 the i = 1. Here the 
early part of the decay of the initial pulse is dominated by 
quasiperiodic peaks with a period approximately equal to 
twice the light travel time between the maximum of the 
exterior effective potential and the origin. Recall that for 
increasingly large compactions, the origin is pushed to 
larger negative values of , thereby increasing the width 
of the potential well. For sufficiently large compactions, 
the wavelength of the incoming pulse is smaller than the 
width of the potential well. Instead of exciting a quasi- 
normal mode, it oscillates between the origin, where it 
is completely reflected, and the maximum of the exter- 
nal Schwazschild potential, where it is partially reflected. 
Every time the pulse reaches the latter, a small part of 
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FIG. 5. Late-time behaviour for the I = and I = 1 modes with C = 0.882. Both modes are in the strongly relativistic 
regime, where the early part of the decay is dominated by almost periodic peaks. In the case of the £ = mode we can also 
see how the exponential decay of quasi-normal modes is recovered in later times, to eventually give rise to the power-law tail. 



it will be transmitted as one of the quasi-periodic peaks 
which we observe at large radii. 

During its travel in the non-zero potential background 
within the well and with each reflection at the barrier, 
the pulse is distorted and becomes broader with time. 
This effect makes the observed pattern quasi-periodic as 
opposed to periodic. Once the typical wavelength be- 
comes comparable with the width of the well, we can 
no longer cleanly separate individual peaks. Simultane- 
ously, the exponential decay found in the intermediate 
regime is recovered, suggesting that once the pulse "fills 
out" the well, the dynamics can again be described by a 
exponentially decaying quasi-normal mode. Ultimately, 
a power-law tail again emerges for very late times. We 
show this behavior for C = 0.882 in Fig. 5. For £ = 1, the 
reappearence of the power-law occurs only once the sig- 
nal has decayed below computer accuracy, which is why 
we are unable to observe this transition. 

An extreme case of the decay pattern described above 
is exhibited in the case of extremely large compactions. 
In this case, as C approaches 8/9, the light travel time 
from the maximum of the exterior potential to the ori- 
gin grows without limit. Thus, the time between the 



initial peaks also grows without bound. As a result, we 
can observe a power- law decay between individual peaks. 
When the next peak arrives, it dominates over the power 
law tail until it decays enough so that the power-law 
once again emerges. However, since the power-law it- 
self drops to increasingly small amplitudes, and since the 
peaks broaden and eventually overlap, the power-law tail 
is temporarily lost, to re-emerge only at very late times, 
when periodic peaks and subsequent exponential decay 
have decayed sufficiently so that the power law tail be- 
comes once again the dominant pulse remnant. The early 
part of this behaviour for I = and for compactions 
C = 0.888 and C = 0.8888 is shown in Fig. 6. 

In the limit C — > 8/9, the light travel time between po- 
tential maximum and origin becomes infinite and hence 
the first pulse will arrive infinitely late. Thus, in the 
Buchdahl limit we recover the characteristic power-law 
decay of a black hole spacetime. 
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FIG. 6. Late-time behavior for 1 = for stars with ex- 
tremely large compactions. The decaying field has enough 
time to enter the power-law tail before the first of the 
quasiperiodic pulses arrives. Note that as the compaction 
increases, the first pulse arrives at increasingly late times. 
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APPENDIX A: CONSERVATION OF ENERGY 

In this appendix we derive a conserved energy integral 
which can be used as a test of the numerical simulation. 

Since the spacetime ([!]) is static, there exists a Killing 
vector £ = d/dt. Thus, a conserved current can be con- 
structed by contracting £ with the stress energy tensor 
T 



J rp £V 



(Al) 



The conservation of J implies J a a — or, by integrating 
over a 4- volume V, 



( 9 1/2 J a ) ,« 



(Fx = 



(A2) 



Applying Gauss's theorem, eq. ( |A2| ) yields 



We study the backscattering of scalar waves off uni- 
form density neutron stars, varying the compaction of 
the neutron star from the Newtonian limit 2M/R <C 1 to 
the ultrarelativistic Buchdahl limit 2M/R = 8/9, and fo- 
cussing on the excitation of quasi-normal modes and the 
emergence of power-law tails. We find that three distinct 
regimes in the compaction can be identified. 

In the Newtonian regime, the response is dominated 
by a reflection of the initial pulse off the center of the 
star, and subsequent power-law falloff. In an intermedi- 
ate regime, more quasi-normal modes of the neutron star 
excited, which decay exponentially and ultimately drop 
in amplitude below the power-law tail, which dominates 
the late-time falloff. In the highly relativistic regime, in 
which the light travel time between the center of the star 
and the maximum of the curvature potential increases to 
large values, the initial wave package bounces back and 
forth between those two and gives rise to quasi-periodic 
signals. Between these quasi-periodic signals, the field 
decays following a power-law. In the ultrarelativistic 
limit, the light travel time reaches infinity, and the ini- 
tial wave signal will never reach the center of the star. 
The modes of neutron stars can therefore no longer be 
excited, and it is in this sense that the late-time radiative 
decay from from ultrarelativistic neutron stars looses all 
its features and again gives rise to power-law tails. By 
tracking the scattered wave amplitude, a distant observer 
can distinguish whether the central object is a black hole 
or a neutron star, and if the latter, can determine the 
compaction of the star. 



J a d 3 Y, a = 



(A3) 



where SV denotes the closed 3-surface enclosing the 4- 
volume V. The surface element is given by d 3 E Q = 
I /'&\g 1 l' 2 £ a f} 1 sdxP dx 1 dx s , where e Q /3 7< 5 is the Levi-Civita 
alternating symbol. In the following we choose V to be 
the volu me d efined by < r < ro and to < t < ii, so 
that eq. (A3) can be written as 



lJr=0 



J°\g\ 1/2 drd9d(f> 



=0 J</>=0 
■2tt 



J r \g\ 1/2 dtd8d<j) 



=0 J (j>—0 



(A4) 



Here the left hand side represents the difference in energy 
interior to a spherical volume between two instants of 
time to and ti, and the right hand side represents the 
integrated flux through the surface of the volume at r — 
ro during that time interval. 

For a massless scalar wave the stress energy tensor is 



T — 



1 

4?r 



1 



(A5) 



so that the relevant components of the flux J M = Tg 
become 



, -2* 2 



-2A 2 i 2 

e + Z5<P,e 



r 2 sin 2 



(A6) 
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FIG. 7. Energy conservation for I = 0, C = 0.85 and for I = 1, G = 0.8. We plot the total energy E(t) contained within 
r = 50, the integrated flux J(t) across r — 50, and the sum, which has to be equal to E(0), a constant. 



and 



J r 



47T 



(A7) 



In the following we will consider an axisymmetric mode 
(rn = 0), for which all derivatives with respect to <p van- 
ish. 

Inserting the fluxes together with the d eter minant of 
the metric \g\ 1 ^ 2 — e* +A r 2 sin6' into eq. (A4) and per- 
forming the integration over <j) then yields 



2 / A-* 2 I „*-A,„2 



drdfir 
= J dtdiir 2 (p_ r ipfle®~ A 



2 



h 

to 

(A8) 



where we h ave used /i = cos 9. The limits of integration 
in eq. (A8) are from r = to ro, from t = to to t\ and 
from fj, = — 1 to 1. 

We can now define E(t) as the energy inclosed within 
a sphere of radius ro at time t, 



E(» - \ J 



drd/jr 



2 / „A-4> 2 



*-A 2 

e <p 



(A9) 



(with the same limits of integration as in eq. |A8| ) , and 
J(t) as the outgoing flux of energy through the surface of 
the same sphere integrated over the surface of this sphere 
and over the time interval t — to, 



J(t) 



*-A 



Eq. ( |AS| ) can now be rewritten 

E(t) + J(t) = E(t Q ) = const 



(A10) 



(All) 



The expression for the energy and flux integrals can be 
further simplified be decomposing <p into spherical har- 
monics 



ip(t, r, 6, 4>) = tp(t, r)Y im (6, i 



(A12) 



as we do in our dynamical evolution. 

For the spherically sym metr ic case £ = 0, loo = 
l/y/4n, the energy integral (A9) becomes 



2 / A-3> - 2 

r ( e <Pfo 



and the flux integral (A10) 

2n J to 



e*- A <p 2 r ) dr (A13) 



2 ~ ~ 4>-A 



dt 



(A14) 
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For £ = 1 a nd m = we have Y\q — \/3/(47r)/i, and 
eqs. (&9) and ( A10 ) now become 



1 f 



1 ( A-4> ~2 i $-A -2\ j 
r (e Vfl+ e L P,r)0' r 

1 fr„ 



^ / e*+Vdr (A15) 
2tt Jo 



and 



i r* 

J to 



(A16) 



Energy conservation in our dy nami cal code can now 
by checked by verifying that eq. (All) is satisfied at all 
times. In Fig. (7) we show the energy and flux integrals 
for examples with I = and I = 1. The sum of the 
two integrals remains constant with an accuracy which 
increases with increasing grid resolution, implying that 
energy is indeed conserved. 
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